{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a2220d36",
   "metadata": {},
   "source": [
    "# medeiros_STORM_analysis\n",
    "\n",
    "Notebook for STORM analysis for Medeiros, et al., 2023."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22f29251",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data locations\n",
    "\n",
    "dataset_num = 0\n",
    "\n",
    "if dataset_num==0:\n",
    "    dbscan_file = '20220517_L1_NMJ1_A3R67_DBSCAN_v2.xlsx'\n",
    "    dbscan_sheet = 'B1_DBSCAN'\n",
    "    max_proj_file = 'MAX_20220517_CacHalo_JF646_L1_A3R67_001.png'\n",
    "else:\n",
    "    dbscan_file = '20220517_L1_NMJ2_A4R67b_DBSCAN_v2.xlsx'\n",
    "    dbscan_sheet = 'B1_DBSCAN'\n",
    "    max_proj_file = 'MAX_20220517_CacHalo_JF646_L1_A4R67b_005.png'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5487a21",
   "metadata": {},
   "source": [
    "## Initialize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6fc6f3f-d548-4615-93cc-36fede3ffa56",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# CGAL package for point->shape analysis\n",
    "from CGAL.CGAL_Kernel import Point_2\n",
    "from CGAL.CGAL_Alpha_shape_2 import Alpha_shape_2, GENERAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "013cc51b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following line will enable interactive plots within the notebook, and\n",
    "# requires ipympl to be installed. You can comment it out to use default\n",
    "# plotting, or try \"%matplotlib qt\" to pop out interactive plots in new\n",
    "# windows.\n",
    "\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a9b147a-ee1b-4698-9920-ed4c69ce64b8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Utility functions\n",
    "\n",
    "def um_to_px(coords, pxum=6.25):\n",
    "    '''\n",
    "    Convert coordinates from microns in the STORM file to \n",
    "    pixels in the confocal image.\n",
    "    \n",
    "    pxum is pixels per micron, found in the metadata\n",
    "    '''\n",
    "    return coords*pxum - 0.5\n",
    "\n",
    "def get_alpha_shape(coords,alpha):\n",
    "    '''\n",
    "    Return the CGAL object for alpha shape from the point cloud.\n",
    "    Note the object contains shape information for all alphas.\n",
    "    \n",
    "    Coords should be an N by 2 array of (x,y) locations.\n",
    "    \n",
    "    Modified from Mrestani, et al 2021\n",
    "    '''\n",
    "    points=[]\n",
    "    for p in coords:\n",
    "        points.append(Point_2(*p))\n",
    "\n",
    "    als = Alpha_shape_2(points, alpha, GENERAL)\n",
    "    als.make_alpha_shape(points)\n",
    "\n",
    "    return als\n",
    "\n",
    "def get_az_area(coords,alpha):\n",
    "    '''\n",
    "    Find the active zone area, defined by a 2D alpha shape\n",
    "    fit to the point cloud. Coords are the point coordinates\n",
    "    in microns, alpha is a float determining the complexity \n",
    "    of the boundary (smaller is more complex).\n",
    "    \n",
    "    Modified from Mrestani, et al 2021\n",
    "    '''\n",
    "    alpha_shape=get_alpha_shape(coords,alpha=alpha)\n",
    "    alpha_shape.set_alpha(alpha)\n",
    "    alpha_shape.set_mode(1)\n",
    "\n",
    "    triangles=[]\n",
    "    for i,h in enumerate(alpha_shape.finite_faces()):\n",
    "        if i==0:\n",
    "            done = h\n",
    "        elif h==done:\n",
    "            break\n",
    "\n",
    "        if not alpha_shape.classify(h):\n",
    "            continue\n",
    "        cur_triangle=[]\n",
    "        for idx in range(3):\n",
    "            p = h.vertex(idx).point()\n",
    "            cur_triangle.append( [p.x(),p.y()] )\n",
    "        triangles.append(cur_triangle)\n",
    "\n",
    "    if len(triangles)>0:\n",
    "        ntriangles=np.ones((len(triangles),3,3))\n",
    "        ntriangles[:,:,1:]=triangles\n",
    "        A=(sum(np.linalg.det(ntriangles))/2)\n",
    "    else: \n",
    "        A=0\n",
    "    return A\n",
    "\n",
    "def plot_img_coords_alpha(Im, alpha_param, cluster_shapes, pts_all):\n",
    "    '''\n",
    "    Plot segmented alpha_shapes and points on top of low-res image\n",
    "    '''\n",
    "    plt.figure()\n",
    "    plt.imshow(Im[:,:,0], cmap = \"gray\")\n",
    "\n",
    "    for cur_shape in cluster_shapes:\n",
    "        cur_shape.set_alpha(alpha_param)\n",
    "        cluster_color = np.random.rand(1,3)\n",
    "        for cur_edge in cur_shape.alpha_shape_edges():\n",
    "            cur_segment = cur_shape.segment(cur_edge)\n",
    "            #print( cur_segment, cur_shape.classify(cur_edge) )\n",
    "            X_from = cur_segment.source().x()\n",
    "            X_to = cur_segment.target().x()\n",
    "            Y_from = cur_segment.source().y()\n",
    "            Y_to = cur_segment.target().y()\n",
    "        \n",
    "            X_px = um_to_px( np.array([X_from,X_to]) )\n",
    "            Y_px = um_to_px( np.array([Y_from,Y_to]) )\n",
    "            plt.plot(X_px,Y_px,color=cluster_color)\n",
    "            plt.scatter(pts_all[:,0],pts_all[:,1],5,color=cluster_color)\n",
    "    plt.axis('image')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b138354a-8c5b-45a9-9229-74e156d4399a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Primary analysis function\n",
    "\n",
    "def get_cluster_stats(dbscan_file, dbscan_sheet, alpha_param, um_to_px):\n",
    "    \n",
    "    # TODO: Remove hard coded columns\n",
    "    data = pd.read_excel(dbscan_file, sheet_name=dbscan_sheet, usecols=['X Position [µm]', 'Y Position [µm]', 'Lateral Localization Accuracy', 'Clusters - ID'])\n",
    "    \n",
    "    #Remove lateral localization accuracy points below 50nm (0.05 microns)\n",
    "    print(f'Number of poorly localized points: {np.sum( data[\"Lateral Localization Accuracy\"] > 0.05)}')\n",
    "    data_loc_acc = data.loc[data[\"Lateral Localization Accuracy\"] <= 0.05]           \n",
    "    del(data)  # Failsafe, TODO: replace with in place operations\n",
    "    \n",
    "    #Sort data by Cluster ID\n",
    "    #data_sorted = data_loc_acc.sort_values(by = ['Clusters - ID'], ascending = True )\n",
    "     \n",
    "    #Pull out datapoints for cluster IDs and calibrate px/um                  \n",
    "    # TODO: Replace in plotting cell\n",
    "    pts_all = um_to_px(data_loc_acc[data_loc_acc['Clusters - ID']>=0].iloc[:,0:3].to_numpy() ) # pixels\n",
    "    pts_all_ids = data_loc_acc[data_loc_acc['Clusters - ID']>=0].iloc[:,3].to_numpy()  #to color points by IDs \n",
    "    \n",
    "    cluster_ids = data_loc_acc.iloc[:,3].unique()  # Check ID column\n",
    "    #cluster_ids = cluster_ids[cluster_ids>=0]  # Reject noise cluster of -1 Column ID if you have it\n",
    "\n",
    "    cluster_stats = np.full( (len(cluster_ids),5), np.nan )\n",
    "    cluster_shapes = []\n",
    "    for kk, cid in enumerate(cluster_ids):\n",
    "        print(f'Working on cluster {cid}')\n",
    "        if cid>=0:\n",
    "            coords = data_loc_acc[data_loc_acc['Clusters - ID']==cid].iloc[:,0:2].to_numpy()\n",
    "            num_coords = len(coords)\n",
    "            cluster_stats[kk,0] = cid  # Identifier\n",
    "            cluster_stats[kk,1] = get_az_area(coords, alpha_param)  # Area, um^2\n",
    "            cluster_stats[kk,2] = 2*np.sqrt(cluster_stats[kk,1]/np.pi)  # Equivalent diameter, um^2\n",
    "            cluster_stats[kk,3] = num_coords / cluster_stats[kk,1]  # Density, num/um^2\n",
    "            cluster_stats[kk,4] = num_coords   # Number of locs in cluster\n",
    "\n",
    "            cluster_shapes.append( get_alpha_shape(coords, alpha_param) ) # Full shape object  \n",
    "     \n",
    "    #cluster_shapes[1].set_alpha(alpha_param)\n",
    "    cluster_stats_pd = pd.DataFrame(cluster_stats,\n",
    "                                    columns = ['Cluster_ID','ClusterArea ( $\\mu m^2$ )','ClusterDiameter ( $\\mu m$ )','Density','# of Locs'])\n",
    "        \n",
    "    return (cluster_stats_pd, cluster_shapes, pts_all)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2adb3fa",
   "metadata": {},
   "source": [
    "## Load and analyze data\n",
    "\n",
    "Expects an Excel workbook in the current directory, with a sheet containing STORM points.\n",
    "\n",
    "Uses CGAL package for alpha shapes to find piecewise linear outlines of each active zone, and compute geometric properties like diameter and density.\n",
    "\n",
    "As `alpha_param` goes to infinity, the alpha shape becomes a convex hull. As `alpha_param` goes to zero, the active zone will be broken into distinct subclusters, eventually with each point its own cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "295a4bbf-1cb5-40f4-86ff-69d4f7edbd76",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_param = 0.015  # Shape complexity\n",
    "\n",
    "cluster_data = get_cluster_stats(dbscan_file, dbscan_sheet, alpha_param, um_to_px)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33e0122c",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_data[0].sort_values(by = ['Cluster_ID'], ascending = True )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59578098",
   "metadata": {},
   "source": [
    "## Visualize\n",
    "\n",
    "Plot the confocal max-projection image, the individual STORM points, and the alpha shapes together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "534f0871-a040-4573-baf3-92c4dc4db104",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load max projection image from raw nd2 video wih Fiji/ImageJ to maintain alignment\n",
    "\n",
    "Im = plt.imread(max_proj_file)\n",
    "\n",
    "plot_img_coords_alpha(Im, 0.0015, cluster_data[1], cluster_data[2]) \n",
    "\n",
    "plt.savefig('Fig1_v1.svg', format='svg', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b60392fd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  },
  "toc-autonumbering": true,
  "toc-showmarkdowntxt": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
